home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 1 / Cream of the Crop 1.iso / PROGRAM / MATRIX04.ARJ / MATDURBN.C < prev    next >
Text File  |  1992-05-22  |  3KB  |  133 lines
  1. /*
  2. *-----------------------------------------------------------------------------
  3. *    file:    matdurbn.c
  4. *    desc:    Levinson-Durbin algorithm
  5. *    by:    ko shu pui, patrick
  6. *    date:
  7. *    revi:    21 may 92 v0.3
  8. *    ref:
  9. *
  10. *       [1] "Fundementals of Speech Signal Processing," Shuzo Saito,
  11. *       Kazuo Nakata, Academic Press, New York, 1985.
  12. *
  13. *-----------------------------------------------------------------------------
  14. */
  15.  
  16. #include <stdio.h>
  17.  
  18. #include "matrix.h"
  19.  
  20. /*
  21. *-----------------------------------------------------------------------------
  22. *    funct:    mat_durbin
  23. *    desct:    Levinson-Durbin algorithm
  24. *
  25. *        This function solve the linear eqns Ax = B:
  26. *
  27. *        |  v0   v1   v2  .. vn-1 | |  a1   |    |  v1   |
  28. *        |  v1   v0   v1  .. vn-2 | |  a2   |    |  v2   |
  29. *        |  v2   v1   v0  .. vn-3 | |  a3   |  = |  ..   |
  30. *        |  ...                   | |  ..   |    |  ..   |
  31. *        |  vn-1 vn-2 ..  .. v0   | |  an   |    |  vn   |
  32. *
  33. *        where A is a symmetric Toeplitz matrix and B
  34. *        in the above format (related to A)
  35. *
  36. *    given:    R = autocorrelated matrix (v0, v1, ... vn) (dim (n+1) x 1)
  37. *    retrn:    x (of Ax = B)
  38. *-----------------------------------------------------------------------------
  39. */
  40. MATRIX mat_durbin( R )
  41. MATRIX R;
  42. {
  43.     int    i, i1, j, ji, p, n;
  44.     MATRIX    W, E, K, A, X;
  45.  
  46.     p = MatRow(R) - 1;
  47.     W = mat_creat( p+2, 1, UNDEFINED );
  48.     E = mat_creat( p+2, 1, UNDEFINED );
  49.     K = mat_creat( p+2, 1, UNDEFINED );
  50.     A = mat_creat( p+2, p+2, UNDEFINED );
  51.  
  52.     W[0][0] = R[1][0];
  53.     E[0][0] = R[0][0];
  54.  
  55.     for (i=1; i<=p; i++)
  56.         {
  57.         K[i][0] = W[i-1][0] / E[i-1][0];
  58.         E[i][0] = E[i-1][0] * (1.0 - K[i][0] * K[i][0]);
  59.  
  60.         A[i][i] = -K[i][0];
  61.  
  62.         i1 = i-1;
  63.         if (i1 >= 1)
  64.             {
  65.             for (j=1; j<=i1; j++)
  66.                 {
  67.                 ji = i - j;
  68.                 A[j][i] = A[j][i1] - K[i][0] * A[ji][i1];
  69.                 }
  70.             }
  71.  
  72.         if (i != p)
  73.             {
  74.             W[i][0] = R[i+1][0];
  75.             for (j=1; j<=i; j++)
  76.                 W[i][0] += A[j][i] * R[i-j+1][0];
  77.             }
  78.         }
  79.  
  80.     X = mat_creat( p, 1, UNDEFINED );
  81.     for (i=0; i<p; i++)
  82.         {
  83.         X[i][0] = -A[i+1][p];
  84.         }
  85.  
  86.     mat_free( A );
  87.     mat_free( W );
  88.     mat_free( K );
  89.     mat_free( E );
  90.     return (X);
  91. }
  92.  
  93. /*
  94. *-----------------------------------------------------------------------------
  95. *    funct:    mat_lsolve_durbin
  96. *    desct:    Solve simultaneous linear eqns using
  97. *        Levinson-Durbin algorithm
  98. *
  99. *        This function solve the linear eqns Ax = B:
  100. *
  101. *        |  v0   v1   v2  .. vn-1 | |  a1   |    |  v1   |
  102. *        |  v1   v0   v1  .. vn-2 | |  a2   |    |  v2   |
  103. *        |  v2   v1   v0  .. vn-3 | |  a3   |  = |  ..   |
  104. *        |  ...                   | |  ..   |    |  ..   |
  105. *        |  vn-1 vn-2 ..  .. v0   | |  an   |    |  vn   |
  106. *
  107. *    domain:    where A is a symmetric Toeplitz matrix and B
  108. *        in the above format (related to A)
  109. *
  110. *    given:    A, B
  111. *    retrn:    x (of Ax = B)
  112. *
  113. *-----------------------------------------------------------------------------
  114. */
  115. MATRIX mat_lsolve_durbin( A, B )
  116. MATRIX A, B;
  117. {
  118.     MATRIX    R, X;
  119.     int    i, n;
  120.  
  121.     n = MatRow(A);
  122.     R = mat_creat(n+1, 1, UNDEFINED);
  123.     for (i=0; i<n; i++)
  124.         {
  125.         R[i][0] = A[i][0];
  126.         }
  127.     R[n][0] = B[n-1][0];
  128.  
  129.     X = mat_durbin( R );
  130.     mat_free( R );
  131.     return (X);
  132. }
  133.